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We report the first large-scale statistical study of very high-lying eigenmodes (quantum states) of 
the mushroom billiard proposed by L. Bunimovich in this journal, 11, 802 (2001). The phase space 
of this mixed system is unusual in that it has a single regular region and a single chaotic region, and 
no KAM hierarchy. We verify Percival's conjecture to high accuracy (1.7%). We propose a model 
for dynamical tunneling and show that it predicts well the chaotic components of predominantly- 
regular modes. Our model explains our observed density of such superpositions dying as E^ 1 ^ 3 (E 
is the eigenvalue). We compare eigenvalue spacing distributions against Random Matrix Theory 
expectations, using 16000 odd modes (an order of magnitude more than any existing study). We 
outline new variants of mesh-free boundary collocation methods which enable us to achieve high 
accuracy and such high mode numbers orders of magnitude faster than with competing methods. 



Quantum chaos is the study of the quantum 
(wave) properties of Hamiltonian systems whose 
classical (ray) dynamics is chaotic. Billiards are 
some of the simplest and most studied exam- 
ples; physically their wave analogs are vibrating 
membranes, quantum, electromagnetic, or acous- 
tic cavities. They continue to provide a wealth of 
theoretical challenges. In particular 'mixed' sys- 
tems, where ray phase space has both regular and 
chaotic regions (the generic case), are difficult to 
analyse. Six years ago Bunimovich described [1] 
a mushroom billiard with simple mixed dynamics 
free of the usual island hierarchies of Kolmogorov- 
Arnold-Moser (KAM). He concluded by antici- 
pating the growth of "quantum mushrooms": it is 
this gardening task that we achieve here, by de- 
veloping advanced numerical methods to collect 
an unprecedented large number n of eigenmodes 
(much higher than competing numerics [2] or mi- 
crowave studies [3]). Since uncertainties scale 
as ra -1 ' 2 , a large n is vital for accurate spectral 
statistics and for studying the semiclassical (high 
eigenvalue) limit. We address three main issues, 
i) The conjecture of Percival [4] that semiclassi- 
cally modes live exclusively in invariant (regular 
or chaotic) regions, and occur in proportion to the 
phase space volumes, ii) The mechanism for dy- 
namical tunneling, or quantum coupling between 
classically-isolated phase space regions, iii) The 
distribution of spacings of nearest-neighbor eigen- 
values, about which recent questions have been 
raised [3]. We show many pictures of modes, in- 
cluding the boundary phase space (the so-called 
Husimi function). 
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FIG. 1: a) Mushroom billiard Q used in this work. The dot- 
ted line shows the reflection symmetry, b) Desymmetrized 
half-mushroom Q' used for mode calculation, and polar coor- 
dinates. Dashed lines meeting at this corner are zeros enforced 
by basis functions. The remaining part of dQ,' is V, comprising 
two pieces: Dirichlet boundary conditions on the parts shown 
as solid, while boundary conditions vary (see text) on the 
dash-dotted vertical line T s . Boundary coordinate q £ [0, L] 
parametrizes F. 



I. INTRODUCTION 

The nature of eigenfunctions of linear partial differ- 
ential operators in the short wavelength, or semiclassi- 
cal, limit remains a key open problem which continues 
to engage mathematicians and physicists alike. When 
the operator is the quantization of a classical Hamilto- 
nian dynamical system, the behavior of eigenfunctions 
depends on the class of dynamics. In particular, hy- 
perbolic dynamics (exponential sensitivity to initial con- 
ditions, or chaos) leads to irregular eigenfunctions, the 
study of which forms the heart of a field known as 'quan- 
tum chaos' [|[ or 'quantum ergodicity' [1, 0|- The planar 
billiard, or particle undergoing clastic reflection in a cav- 
ity C M 2 , is one of the simplest examples. Billiards 
exhibit a menagerie of dynamical classes [1, [l(| rang- 
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ing from complete integrability (ellipses and rectangles) 
to complete ergodicity (e.g. Sinai billiard). Bunimovich 
introduced the 'mushroom' billiard [H, [ll| with the nov- 
elty of a well-understood divided phase-space comprising 
a single integrable (KAM) region and a single ergodic re- 
gion [561 ] - As seen in Fig.[T^i, the mushroom is the union 
of a half-disk (the 'hat') and a rectangle (the 'foot'); only 
trajectories reaching the foot are chaotic. The simplic- 
ity of its phase space has allowed analysis of phenomena 
such as 'stickiness' (power-law decay of correlations) in 
the ergodic region [l2|, [Hj] . 

The quantum-mechanical analog of billiards is the 
spectral problem of the Laplacian in Q with homoge- 
neous boundary conditions (BCs). Choosing Dirichlct 
BCs (and units such that h = 2m = 1) we have 

- A(f>j = Ej(f>j in fl, (1) 
4> 3 = ondtt. (2) 

This 'drum problem' has a wealth of applications 
throughout physics and engineering [14| . Eigenfunctions 
(or eigenmodes, modes) (f>j may be chosen real- valued and 
orthonormalized, ((f>i,<j)j) := ^ n 4>i(r)<j)j(r)dr — Sij, 
where dr :— dxdy is the usual area element. 'Energy' 
eigenvalues E\ < E2 < E3 < ■ ■ ■ — ► 00 may be written 
Ej = k 2 , where the (eigen)wavenumber kj is 2ir divided 
by the wavelength. 

Traditional numerical methods to compute eigenval- 
ues and modes employ finite differences or finite elements 
(FEM). They handle geometric complexity well but have 
two major flaws: i) it is very cumbersome to achieve high 
convergence rates and high accuracy and ii) since sev- 
eral nodes are needed per wavelength they scale poorly 
as the eigenvalue E grows, requiring of order E degrees 
of freedom (e.g. for the mushroom deMenezes et. al. [l5j 
appear limited to j < 400). The numerical difficulty is 
highlighted by the fact that analog computation using 
microwave cavities is still popular in awkward geome- 
tries nm. 

In contrast we use boundary-based methods, as ex- 
plained in Sec. [TTJ These i) achieve spectral accuracy, 
allowing eigenvalue computations approaching machine 
precision as exhibited for low- lying modes in Sec. lIIIl and 
ii) require only of order E 1 ! 2 degrees of freedom (with 
prefactor smaller than boundary integral methods [17|). 
Furthermore at high E we use an accelerated variant, 
the scaling method QjJ [2(| , which results in another 
factor of order E 1 / 2 in efficiency. These improvements 
allow us to find large numbers of modes up to j ~ 10 5 ; in 
Sec. lIVI we show such modes along with their Husimi (mi- 
crolocal) representations on the boundary. Visualization 
of modes can be an important tool, e.g. in the discovery 
of scars (2lj . 

We are motivated by a growing interest in quantum er- 
godicity [12, [Hj|. For purely ergodic billiards, the Quan- 
tum Ergodicity Theorem [|l [Si [H, [27[ (QET) states 
that in the E — > 00 limit almost all modes become 
equidistributed (in coordinate space, and on the bound- 
ary phase space [1^, However no such theorem ex- 




FIG. 2: The tension t m (E) plotted as a function of trial 
eigenvalue parameter E, for the half-mushroom with Dirichlet 
boundary conditions. The minima indicate the eigenvalues of 
this domain. Close to E — 44 there is a cluster of two eigen- 
values. 



ists for mixed billiards, thus numerical studies are vital. 
It is a long-standing conjecture of Percival [H that for 
mixed systems, modes tend to localize to one or another 
invariant region of phase space, with occurence in pro- 
portion to the phase space volumes, and that those in er- 
godic regions are equidistributed. (This has been tested 
in a smooth billiard [30] , and recently proved for certain 
piecewise linear quantum maps [3l[). We test the conjec- 
ture via a matrix element (|I0p sensitive to the boundary 
(for numerical efficiency); we then can categorize (almost 
all) modes as regular or ergodic. We address two issues 
which have also been raised by recent microwave exper- 
iments in the mushroom [3J. i) The mechanism for dy- 
namical tunneling [32| is unknown (although it has been 
studied in KAM mixed billiards [3a]). In Sec. ED we pro- 
pose and test a simple model for coupling strength (re- 
lated to [13]) which predicts observed features of matrix 
element distributions, ii) The level-spacing distribution, 
conjectured to be a universal feature 0, [3||, is studied 
in Sec. I VII where we also examine spacing distributions 
for regular and ergodic subsets of modes. Note that we 
use an order of magnitude more modes than any exist- 
ing experiment or study. Finally we draw conclusions in 
Sec. En] 



II. NUMERICAL METHODS 

In this section we outline the numerical methods that 
make our investigation possible; the reader purely inter- 
ested in results may skip to Sec. Mil 



A. The Method of Particular Solutions 

Our set of basis functions, or particular solutions, 
{C"( r )}n=i - A r satisfy — A£„ = E£ n at some trial eigen- 
value parameter E, but do not individually satisfy 
The goal is now to find values of E such that there exists 
nontrivial linear combinations xi^i + 2:2^2 + ■ ■ ■ + xn(,n, 
which are small on the boundary. These are then hope- 
fully good approximations for an eigenfunction. 
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Let us make this precise. We define the space Ti-(E) of 
trial functions at a given parameter E as 

H{E) =Span{a,.--,6r}- 

If we denote by ||u||aa and the standard L 2 -norm 

of a trial function u £ Ti-(E) on the boundary <9f2 and 
in the interior f2 we can define the normalized boundary 
error (also called the tension) as 



t[u] := 



Mian 

Nln 



(3) 



It is immediately clear that t[u] =0 for u £ H-(E) if and 
only if u is an eigenfunction and E the corresponding 
eigenvalue on the domain f2. However, in practice we 
will rarely achieve exactly t[u] = 0. We therefore define 
the smallest achievable error as t m (E) :— min ue -^(B) t[u]. 
This value gives us directly a measure for the error of 
an eigenvalue approximation E, namely there exists an 
eigenvalue Ej such that 



E, 



< Ct m (E), 



(4) 



where C is an O(l) constant that only depends on the 
domain f2. This result is a consequence of error bounds 
of Moler and Payne [36L l37j| . Hence, by searching in E for 
minima of t m (E) we find approximate eigenvalues with 
relative error given by a constant times t m (E). Fig. [5] 
shows such a plot of t m (E) for our mushroom domain. 

The implementation of this Method of Particular So- 
lutions (MPS) depends on i) basis set choice, and ii) how 
to evaluate t m (E). The former we address in the next 
section. The latter requires a set of quadrature points 
{yi}i=i - M on which to approximate the boundary inte- 
gral ||u||an. One must take into account that Helmholtz 
basis sets tend to be ill-conditioned, that is, the M x N 
matrix A with entries Ai n := £ n (yi) becomes numerically 
rank-deficient for desirable choices of N. The tension 
t m (E) can then be given by the square-root of the lowest 
generalized eigenvalue of the matrix pair (A T A, B T B), 
or by the lowest generalized singular value of the pair 
(A,B), where B is identical to A except with the re- 
placement of {yi} by interior points fl4l. Il9l . l38j . These 
different approaches are discussed in [39j|. Here, we use 
the generalized singular value implementation from [39| . 
which is highly accurate and numerically stable. We note 
that these methods are related to, but improve upon, the 
plane wave method of Heller • 



B. Choice of basis functions 

In order to obtain accurate eigenvalue and eigenfunc- 
tion approximations from the MPS it is necessary to 
choose the right set of basis functions. In this section 
we propose a basis set that leads to exponential conver- 
gence, i.e. errors which scale as e~ cN for some c > 0, as 
N the number of basis functions grows. 



To achieve this rate we first desymmetrize the problem. 
The mushroom shape Q is symmetric about a straight 
line going vertically through the center of the domain 
(see Fig. [T]). All eigenmodes are either odd or even sym- 
metric with respect to this axis. Hence, it is sufficient 
to consider only the right half, f2'. The odd modes are 
obtained by imposing zero Dirichlet boundary conditions 
everywhere on the boundary dfl' of the half mushroom. 
The even modes are obtained by imposing zero Neumann 
conditions on the symmetry axis T s and zero Dirichlet 
conditions on the rest of dfl'. 

Eigenfunctions of the Laplacian are analytic every- 
where inside a domain except possibly at the bound- 
ary [4l|. Eigenfunctions can be analytically extended 
by reflection at corners whose interior angle is an inte- 
ger fraction of it [lij]. The only singularity appears at 
the reentrant corner with angle 3n/2 (where dashed lines 
meet in Fig. [TJ)). Close to this corner any eigenfunction 
(j>j can be expanded into a convergent series of Fourier- 
Bessel functions of the form 



<M r > 6 = ^2 a kJ^{kjr) sin^6», 



(5) 



71=1 



where the polar coordinates (r, 9) are chosen as in Fig.[T|D. 
The function J a is the Bessel function of the first kind of 
order a. 

The expansion (J3J) suggests that the basis set := 
J^n^kr) sin where k 2 = E, might be a good choice 
since these functions capture the singularity at the reen- 
trant corner and automatically satisfy the zero bound- 
ary conditions on the segments adjacent to this corner 
(dashed lines in Fig. [TJd). Hence, we only need to min- 
imize the error on the remaining boundary T which ex- 
cludes these segments. The boundary coordinate q £ 
[0,L] parametrizes T; its arc length is L = 3(1 + 7r/4). 
This Fourier-Bessel basis originates with Fox, Henrici and 
Moler [i^l for the L-shaped domain; we believe it is new 
in quantum physics. In [43[ the convergence properties 
of this basis set are investigated and it is shown that for 
modes with at most one corner singularity the rate of 
convergence is exponential. Indeed, in practice we find 
t m (Ei) = 0(e~ cN ) for some c > as the number N 
of basis functions grows. Hence, for the minimum E of 
t m (E) in an interval containing E\ it follows from (|3|) 
that 



\E-E X \ 
Ei 



<Ct m {E)<Ct m {E l ) = 0{e~ cN ), 



which shows the exponential convergence of the eigen- 
value approximations E to E\ for growing N. 



C. Scaling method at high eigenvalue 

For all odd modes apart from the lowest few we 
used an accelerated MPS variant, the scaling method 
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FIG. 3: The first 10 odd (a) and even (b) modes of the mush- 
room shape, shown as density plots. Eigenvalue increases 
rightwards from the top left. White corresponds to positive 
and black to negative values. 



18, 19, 20] , using the same basis as above (to our knowl- 
edge the scaling method has not been combined with a re- 
entrant corner-adapted basis before now). Given a center 
wavenumber fco and interval half-width Afc, the scaling 
method finds all modes <pj with kj € [fco — Afc, fco + Afc]. 
This is carried out by solving a single indefinite general- 
ized eigenvalue problem involving a pair of matrices of the 
type A T A discussed above. The 'scaling' requires a choice 
of origin; for technical reasons we are forced to choose the 
singular corner. Approximations to eigenvalues lying in 
the interval are related to the matrix generalized eigen- 
values, and the modes to the eigenvectors. The errors 
grow [lj| as \kj — fco| 3 , thus the interval width is deter- 
mined by the accuracy desired; we used Afc = 0.1 which 
ensured that t m {E) errors associated with the modes 
rarely exceeded 3 x 10~ 4 . Since the search for minima 
required by the MPS has been avoided, and on average 
0{k) modes live in each interval, efficiency per mode is 
thus 0(k) = 0(£ 1/2 ) greater than the MPS. By choosing 
a sequence of center wavenumbers fco separated by 2Afc, 
all modes in a large interval may be computed. Rather 
than determine the basis size N by a convergence crite- 
rion as in Sec. Ill Bl for E > 10 3 we use the Bessel func- 
tion asymptotics: for large order J a (x) becomes expo- 
nentially small for x/a < 1 (the turning point is x = a). 
Equating the largest argument kR (with R = 3/2) with 
the largest order 2iV/3 gives our semiclassical basis size 
N w 9fc/4 = 0{E X / 2 ). 

We are confident that the scaling method finds all odd 



modes in a desired eigenvalue window. For instance we 
compute all 16061 odd symmetry modes with kj < 300, 
using 1500 applications of the scaling method (at fco = 
0.1, 0.3, 299.9). This computation takes roughly 2 
days of CPU time [57[ ■ We verify in Fig. [5] that there 
is zero mean fluctuation in the difference between the 
(odd) level-counting function N(k) := #{j : kj < k} and 
the first two terms of Weyl's law @ , 

where \dfl'\ is the full perimeter of the half mushroom do- 
main. Note that there is no known variant of the scaling 
method that can handle Neumann or mixed BCs, hence 
we are restricted to odd modes. It is interesting that 
the method is still not completely understood from the 
numerical analysis standpoint [lg, [li| [2(| • 

In applying the scaling method to the mushroom, 
the vast majority of computation time involves evalu- 
ating Bessel functions J a (x) for large non-integral a and 
large x. This is especially true for producing 2D spa- 
tial plots of modes as in Fig. for which of order 10 9 
evaluations are needed (1 hr CPU time). We currently 
use independent calls to the GSL library [44[ for each 
J a (x) evaluation. This is quite slow, taking between 
0.2 and 50 (is per call, with the slowest being in the 
region a < 50, 10 2 < x < 10 3 . However, we note 
that Steed's method [45|, |46j], which is what GSL uses 
in this slow region, is especially fast at evaluating se- 
quences J a (x), J a -i(x), J a -2(x), . . ., and that since a is 
a multiple of a rational with denominator 3, only 3 such 
sequences would be needed to evaluate all basis functions 
{£m( r )}m=i - Af at a given location r. We anticipate at 
least an order of magnitude speed gain could be achieved 
this way. 




FIG. 4: The 10 odd modes of the mushroom whose eigen- 
wavenumbers lie in the range 90 < kj < 90.35, at mode num- 
ber about j » 1430. Intensity \<j>j\ 2 is shown with zero white 
and larger values darker. 
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11.50790898 
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5.497868889 


2 


25.55015254 


2 


13.36396253 


3 


29.12467610 


3 


18.06778679 


4 


43.85698300 


4 


20.80579368 


5 


44.20899253 


5 


32.58992604 


6 


53.05259777 


6 


34.19488964 


7 


55.20011630 


7 


41.91198264 


8 


66.42332921 


8 


47.37567140 


9 


69.22576822 


9 


54.62497098 


10 


82.01093712 


10 


65.18713235 



TABLE I: Tables of a) lowest 10 odd and b) lowest 10 even 
eigenvalues of the mushroom. All digits shown are believed 
to be correct. 



III. LOW EIGENVALUE MODES 

In this section we present highly-accurate results for 
the first few even and odd modes. Odd modes are 
obtained by solving the eigenvalue problem with zero 
Dirichlet boundary conditions on the half mushroom 
from Fig.[]J>, using the MPS, by locating minima in the 
tension function of Fig. [2l In Tabled the eigenvalues are 
listed to at least 10 significant digits, and in Fig. the 
corresponding modes are plotted. We emphasize that it 
is the exponential convergence of our basis that makes 
such high accuracies a simple task. 

For even modes we impose Neumann BCs on T s and 
Dirichlet BCs on the remaining part of T. This was 
achieved in the MPS by modifying the tension function 
([3]) to read 




where the normal derivative operator on the boundary is 
d n :— n ■ V, the unit normal vector being n. Table Up, 
gives the smallest 10 even modes on the mushroom bil- 
liard, and the corresponding modes are plotted in Fig.(3jD. 

Although we are far below the semiclassical regime we 
already see properties of the underlying classical dynam- 
ical system. For example, the 8 th odd and the 6 th even 
mode live along a caustic and therefore show features of 
the classically integrable phase space while the 7 th odd 
and 10 th even mode already shows features of the clas- 
sically ergodic phase space. For comparison, in Fig. [4] 
we show some odd modes with intermediate eigenvalues 
of order 10 4 (odd mode number of order 10 3 ), a similar 
quantum number to that measured in a microwave cav- 
ity by Dietz et. al. [|[. As these authors noted, modes at 
this energy usually live in either the integrable or to the 
ergodic regions of phase space; we pursue this in detail 
in Sec. El 
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FIG. 5: Difference between the mode counting function N(k) 
and the two-term Weyl's prediction JVweyi(fc) defined by ||B). 
for the 16061 odd-symmetric eigenvalues with kj < 300. 
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FIG. 6: Poincare surface of section (PSOS), that is, the clas- 
sical phase space in boundary coordinates q (as shown in 
Fig- Eb) an d p (sin of incidence angle). Vertical dashed lines 
shows location of corners. The vertical dotted line shows lo- 
cation of q c the focal point of the hat. The dark line shows 
the border of integrable phase space; note that q = 2 corre- 
sponds to the smallest possible caustic for integrable phase 
space. Families of orbits defined by constant angular momen- 
tum are shown by lines in the integrable region. Note that 
they exchange vertical ordering at the corner, as indicated by 
their grayscale color labeling. 



IV. BOUNDARY AND HUSIMI FUNCTIONS 

We choose a Poincare surface of section (PSOS) @ de- 
fined by Birkhoff coordinates (q,p) G T x [—1,1] =: Z, 
where q is the boundary location as before (see Fig. [TJd) 
and p the tangential velocity component, in the clock- 
wise sense, for a unit speed particle. (If the incident 
angle from the normal is 8 then p = sin 6). The struc- 
ture of this PSOS phase space is shown in Fig. [51 Our 
choice (which differs from that of Porter et. al. [ll]]) is 
numerically convenient since it involves only the part of 
the boundary on which matching is done (Sec. [IT]). De- 
spite the fact that it does not cover the whole boundary 
<9£1', it is a valid PSOS since all trajectories must hit T 
within bounded time. 

Integrable phase space consists of precisely the orbits 
which, for all time, remain in the hat [l[ but which never 
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FIG. 7: Intensity of boundary normal-derivative functions 
\dn4>j(q)/kj\ 2 , plotted vs boundary coordinate g on the hori- 
zontal axis and odd mode number j £ [1, 600] on the vertical. 
The density plot shows white as zero, and larger values darker. 



come within a distance b/2 from the center point q c [58j | . 
Simple geometry shows that the curved boundary be- 
tween ergodic and integrable regions consists of points 
(q,p) satisfying 



b/2 



p 



for p 2 <p 2 == 1-7^2- ( 8 ) 



AR 2 ' 



For our shape, q c — a + b/2 = 3/2, p\ = 8/9. In the 
domain g g [q c + R, L] the boundary occurs at the lines 
p = ±b/2R = ±1/3. Successive bounces that occur on 
r are described by the PSOS billiard map / : Z — > Z. 
Any such Poincare map is symplectic and therefore area- 
preserving ||. 

The quantum boundary functions d n <j)j{q) for q S [0, L] 
are convenient and natural representations of the modes. 
Note that they are not L 2 (dft) normalized; rather they 
are normalized according to a geometrically- weighted L 2 
boundary norm via the Rellich formula (see [2G. 1471]) 



(r-n)\d n h\ 2 dq = 2E 



j > 



(9) 



where r(q) is the location of boundary point q relative to 
an arbitrary fixed origin. Fig. [7] shows the intensities of 
the first 600 odd boundary functions. Features include an 
absence of intensity near the corners (over a region whose 
size scales as the wavelength). The region 3 < q < L, 
in which phase space is predominantly integrable, has a 
more uniform intensity than < q < 2, which is exclu- 
sively ergodic. The region 2 < q < 3 is almost exclusively 
integrable, but is dominated by classical turning-points 
corresponding to caustics; these appear as dark Airy-like 
spots. In 1/2 < q < 3/2 there are horizontal dark streaks 
corresponding to horizontal 'bouncing-ball' (BB) modes 
in the foot. Finally, a series of slanted dark streaks is 
visible for 3/2 < q < 2: these interesting fringes move as 
a function of wavenumber and we postpone analysis to a 
future publication. 

In Fig. Q3] we show a sequence of 20 much higher modes 
with consecutive eigenvalues near wavenumber k = 500 
(eigenvalue E — 2.5 x 10 5 ). These modes are a subset 
of the modes produced via a single generalized matrix 
eigenvalue problem (of size N w 1200) using the scaling 
method at fco = 500. The full set of 77 modes (evalu- 
ating boundary functions) took only 20 mins CPU time. 
Typical tension t m (E) values were below 10 -3 . Naively 
applying ^ we would conclude only about 3 relative 
digits of accuracy on eigenvalues. However, it is possible 
to rigorously improve this bound by factor 0(E^ 2 ) [3^ |. 
giving about 6 digits. 

Fig. [5] shows the 14th in the sequence in more detail. 
The corresponding boundary function is shown in Fig. [5^, 
along with the intensity, and its Husimi distribution. The 
Husimi distribution is a coherent-state projection of the 
mode onto the PSOS phase space (see App. IA"|) , The 
choice of the aspect ratio a is somewhat arbitrary but 
it is expected [21j that phase space structures have spa- 
tial scale 0(fc -1 / 2 ), so we chose a scaling similar to this: 
with k = 500 we used a = 0.076. By comparing to 
the phase space (Fig. [6]) we see localization to the er- 
godic region. The only part of ergodic phase space not 
well covered contains BB modes in the foot (the white 
'box'). A scar is also visible as the 9 darkest spots: 4 
pairs of spots surrounding the white box correspond to 
4 bounces in the foot, and a single spot at q « 5 corre- 
sponds to a normal-incidence bounce off the circular arc. 
By contrast, Fig. [^b shows the boundary function of a 
mode living in the regular region (the 15th in Fig. [8J; 
the energy-shell localization is clear. The full set of 20 
Husimi functions is shown in Fig. 1151 We remind the 
reader that in purely ergodic systems boundary functions 
obey the QET [H, [29| with almost every d n (j)j/kj tend- 
ing to an invariant Husimi density of the form C yl — p 2 . 
We might expect a similar result for the ergodic subset 
of modes in the ergodic phase space of the mushroom. 
However, Fig. [TBI highlights that, despite being at a high 
mode number of roughly 4 x 10 4 , we are still a long way 
from reaching any invariant density: the 7 ergodic modes 
have highly non-uniform distributions. 



FIG. 8: High-energy eigenmode with kj = 499.856 • • ■ , at 
around odd mode number j ~ 45000. This mode appears to 
live in the ergodic region. 



V. PERCIVAL'S CONJECTURE AND 
DYNAMICAL TUNNELING 

In the small set of 20 high-lying modes discussed above, 
Percival's conjecture holds: modes are either regular or 
chaotic but not a mixture of both. We will now study this 
statistically with a much larger set, the first n = 16061 
odd modes corresponding to < kj < 300. Since the 
PSOS phase space in0<g<3/2is ergodic for all p, the 
following 'foot-sensing' quadratic form, or diagonal ma- 
trix clement, is a good indicator of an ergodic component: 

/j := 2^/ 3/2(l " n)|5 "^ |2d<? ' (10) 
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FIG. 9: a) the mode of Fig. [8] Husimi distribution 
Hd n 4, j: a(q,p) defined by (|A4jl (top), density plot of \d„(f>j\ 2 
(middle), and graph of d„<j>j (bottom). Note the q coordinate 
is common to the three plots, b) Similar representation of the 
next highest mode at kj = 499.858, the 15 th in the sequence 
of Fig. 1141 which lives in the regular region. 



where, as Fig. [T}d shows, r • n takes the value 1 for < 
q < 1/2, and 1/2 for 1/2 < q < 3/2. (The weighting by 
r ■ n is chosen to mirror ([9|); scaling by Ej is necessary 
for a well-defined semiclassical limit |28l l48l|). 

The observed distribution of fj is shown in Fig. [TUb . 
The main feature is a cluster around 0(1) (we associate 
with ergodic modes) and a wider distribution of smaller 
values (predominantly regular modes). We have tested 
that the apparent cluster lying roughly from 10~ 14 to 
10~ 9 is merely an artifact reflecting the size of numerical 
errors in d n <t>j'- the key point is that there is a continuum 
of values (see errorbars in Fig. [TUb ) which extends from 
0(1) down to exponentially small values. Roughly 0.75% 
of the total number of modes fall within each decade from 
10~ 2 to 10~ 8 . We believe that in the absence of numerical 
errors a similar distribution would extend down many 
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FIG. 10: a) Histogram of the logarithm of fj, the 'foot- 
sensing' matrix element (|10p . for the first 16061 odd modes. 
Errorbars show counts of fj lying in each decade (errors as- 
suming independent counts), on a vertical scale magnified by 
a factor 15. b) Fraction of modes with 10~ 8 < fj < 10~ 2 lying 
in logarithmically-spaced Ej intervals (errorbars), compared 
to power law _E -1 ' 3 (solid line). 



tens of orders of magnitude. 

Percival's conjecture would imply that the sequence 
{/j}j=i— oo has (for all but a set of vanishing measure) 
two limit points: zero (for regular modes), and some pos- 
itive constant (for ergodic modes). Even though most 
mode numbers are large (~ 10 4 ) the upper cluster still 
has a wide standard deviation of 0.1 (its mean is 0.39); 
this is in line with our recent work confirming the slow al- 
gebraic semiclassical convergence of matrix elements [2(| ■ 

We would like to test whether the relative mode fre- 
quencies of regular vs ergodic modes are in proportion to 
the corresponding classical phase-space volumes. We cat- 
egorize modes by defining them as 'regular' if fj < 0.1. 
This choice of cutoff value is necessarily a compromise 
between lying below the whole ergodic peak yet captur- 
ing the full dynamic range of regular modes. This gives 
a fraction a rog := n re Jn = 7178/16061 = 0.4469 •• • of 
regular modes, which is only 1.7% less than the inte- 
grable phase space fraction jU reg = 0.4549 • • • (computed 
in App. [Bj. Assuming that each regular mode counted 
arose randomly and independently due to some under- 
lying rate (fraction of level density), we may associate 
a standard error of y/n rt . e (n — n rc(! )/n 3 = 0.004 with 
the measured fraction. Thus the discrepancy is only 2 
sigma, not inconsistent with the (null) hypothesis that 
a reg = j« reg . To check whether this result persists semi- 
classically we computed a smaller set of n = 615 high- 
lying modes sampled from the range 500 < kj < 750, up 
to mode number j sa 10 5 , and found a rt , s = 0.441 ±0.015, 
again consistent with Percival's conjecture. 



A. Results and model for dynamical tunneling 

The continuum of matrix element values in Fig. llOb , is a 
manifestation of dynamical tunneling [32j | , quantum cou- 
pling between regular and ergodic invariant phase space 
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FIG. 11: a) Distribution of matrix elements fj for regular 
modes, as function of eigenvalue Ej. b) Corresponding rates 
7j predicted by the model in the text, c) Zoom of a), showing 
fj (dots) connected by vertical lines to best-matching values 
of C7j (circles labeled by disc quantum numbers m, n), using 
constant c = 15. d) ratio fj /jj for all matched pairs. 



regions. This has recently been seen in mushroom mi- 
crowave cavity modes Q, and these authors raised the 
question as to the mechanism for tunneling in this shape. 
We address this by proposing and numerically testing a 
simple such model. First we notice that the density of 
In fj is roughly constant (in the range fh > 10~ 8 where 
numerical errors are negligible). This suggests a cou- 
pling strength which is the exponential of some uniformly- 
distributed quantity. We may ask whether this density 
is dependent on eigenvalue magnitude (energy): Fig. llOb 
shows that the density appears to die as i? -1 / 3 , consis- 
tent with the expectation that all fj values for regular 
modes vanish in the semiclassical limit. 

Our model is to assume that fj values are controlled 
by a matrix element jj giving the rate of dynamical tun- 
neling from the regular to the ergodic region. Each reg- 
ular mode closely approximates an (n, m)-mode of the 
quarter disc, which are the product of angular function 
(2/y^r) shimf? and radial function 



1pmn(r) 



V2 



J m {k^ 



(ii) 



where n = 1,2,... is the radial mode number and 
m = 2,4,... the angular mode number, and k mn R is 
the argument of the n th zero of the J m Bessel function. 
Quarter-disc eigenwavenumbers are k mn . The normal- 
ization is §q \ipmn{f)\ 2 rdr = 1. A wavepacket initially 
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launched from such a disc mode will, in the mushroom, 
leak into the ergodic region due to the openness of the 
connection into the foot. We take the rate proportional 
to the probability mass of ip mn 'colliding' with the foot, 



b/2 

1] ■= I \ip mn (r)\ 2 rdr 
h 



(12) 



J2{m + l + 2l) 



1=0 



Jm+l+2l I 



where we used [49, Eq. 11.3.2] to rewrite the integral. 
This model is similar to that proposed recently by Backer 
et. al. 34] (in our case the 'fictitious integrable system' 
is the quarter-disc), jj is exponentially small only when 
the Bessel function turning point lies at radius greater 
than 6/2; at eigenvalue E this occurs for b\fEj(2m) < 1. 

We compare in Fig. [TTk ) and b) fj values for regu- 
lar against jj values computed using all relevant (m, n) 
quantum numbers for the quarter-disc. It is clear that 
although the densities are similar, fj is irregularly dis- 
tributed whereas jj values fall on a regular lattice. How- 
ever, upon closer examination there is a strong corre- 
lation. We attempted to match each disc mode (m, n) 
seen in panel b) with its corresponding mushroom mode 
j as seen in panel a); in most of the 1051 cases there 
was a very clear match, with relative eigenvalue differ- 
ence \Ej — k^nl/Ej < 10~ 4 in 90% of the cases, and 

\ E 3 ~ k Ln\l E 3 < 3 x 10_6 in 74% of cases - ( Note that ' 
although it is not needed for our study, it would likely 

be possible to improve the fraction matched using data 



from d„ 



As shown in panel c), /,■ values are quite 



correlated with the jj values of their matched mode. 
Note that an overall prefactor of c = 15 was included 
to improve the fit. The resulting ratio fj/^ij is shown 
in panel d), and has a spread of typically a factor 10 2 . 
Since this is much less than the spread of 10 8 in the origi- 
nal matrix elements, this indicates that the above model 
is strongly predictive of dynamical tunneling strength, 
mode for mode. We suggest the remaining variation, 
and the value of c, might be explained by varying eigen- 
value gaps (resonant tunneling) between quarter-disc and 
ergodic modes (such variation is discussed in (33[), al- 
though this is an open question. Also in this simple 
model it is clear from the ^-dependence in panel d) that 
there are algebraic prefactors that should be included in 
a more detailed model. 

Using the model we may predict the decay ~ E^ 1 ^ 3 
in the density of fj values reported above, by returning 
to the sum in (|12p. For regular modes where jj <C 1, 
the Bessel functions in the I > 1 terms have turning 
points successively further away from 6/2, thus the sum 
may be approximated by the / = term (this has been 
checked numerically). We make the approximation that 
the turning point is close to 6/2, that is e <C 1, where 



£ := 1 



bk r , 



2m 



(13) 



drop algebraic prefactors. In (I12[) using Debye's asymp- 
totics for the Bessel function [49j, Eq. 9.3.7] and keeping 
leading terms for small e > gives 
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ln 7j 



(const) + — In e 



i , 4 ^ 3/2 

lnfc ro „ + — me 3 ' 2 . 

(This can be interpreted as the tail of the Airy ap- 
proximation to the Bessel). For fixed e « 1 we need 
keep only the last term as m — > oo. Fixing m while 
increasing n by 1 causes a small wavenumber change 
km,n+i — k mn w tt/(pqR), causing via (fT3"|) a change 
Ae w — irb/ (2mpoR), which in turn causes via (fl"4|) a 
change 



nV2e/(p Q R) 



3g_ 
PoR \2m 



1/3 



(15) 



where in the last step we expressed e in terms of the 
asymptotic for g. Realising that, for e *C 1 we have 
m rs bk mn /2 = by/E/2, and that adjacent curves of 
constant m in the (E, g)-plane are separated in E by 
AE w 8V£/6, gives our result, the density of points in 
the (E, <7)-plane, 



d(E,g) 



1 



PoR ( b 4 



|A 3 A£;| 



8tt \3gE 



1/3 



(16) 



We focus on the exponentially small behavior of jj and 



Recall that Fig. [TTk.b.c illustrate the (E, g)-plane. In 
Fig. 110b small dynamic range and counting statistics pre- 
vents this weak dependence of density on g from being de- 
tected. However the main conclusion from is that the 
density of 7j (and hence fj) values lying in any fixed in- 
terval scales asymptotically as E^ 1 ! 3 , in agreement with 
Fig.[TOb. 



VI. LEVEL SPACING DISTRIBUTION AND 
LEVEL DENSITY FLUCTUATION 

We show the nearest-neighbor spacing distribution 
(NNDS) of the complete set of the first n = 16061 eigen- 
values of odd-symmetric modes Ej with kj < 300, in 
Fig. 112b . Spacings were unfolded in the standard way 
[35j . thus a histogram of Sj := (Ej+i — Ej)/E, where E 
is the mean level spacing, was collected. This is compared 
in the figure against the Berry-Robnik prediction (50j 
for a mixed system with a single regular component (of 
phase-space fraction /i rcg = 0.4549 • • • ) and single ergodic 
component. The agreement is excellent, with deviations 
consistent with the standard error for each bin count. In 
their recent work Dietz et. al. 0] claim that there is a dip 
in the NNDS around s — 0.7 associated with supershell 
structure in the hat (two periodic orbits of close lengths). 
Their choice of mushroom shape differs from ours only 
in the foot. Our results, computed using over 16 times 
their number of levels, show no such dip. This suggests 
that their observed dip is a statistical anomaly, or that 
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a) all levels 



b) regular 




FIG. 12: Nearest-neighbor spacing distributions (NNDS) p(s) 
for the 16061 modes with kj < 300, estimated via a histogram 
with bins of width As = 0.125. Data (in counts per bin) are 
shown by bars. Predictions are shown by dots, with ±1 stan- 
dard error (solid lines above and below), a) all modes, vs 
Berry-Robnik formula, b) regular modes, vs Poissonian for- 
mula e~ s , c) ergodic modes including BB modes, vs Wigner's 

approximate GOE formula -|se _7rs and d) ergodic modes 
with BB modes removed using the categorization in Sec. IVII 
vs the same. 
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FIG. 13: Absolute value of the Fourier transform p(l) of the 
density of states p(k) := YlTLj — kj), vs orbit length I, 
for a) all eigenvalues lying below kj < 300, b) regular modes 
only, c) ergodic modes only. 



it does not carry over to the rectangular-foot mushroom 
and therefore is not associated with the hat. 

In order to study this further we computed the partial 
NNDS associated with regular or ergodic modes, cat- 
egorized using the method of Sec. [V] Regular modes 
(Fig. [T2"b) fit the Poisson level spacing distribution well. 



Ergodic modes (Fig. [P2b l fit Wigner's standard approx- 
imate form for the GOE distribution reasonably well, 
however there are visible deviations: the data system- 
atically favors small spacings s < 0.75 while disfavoring 
intermediate spacings 0.75 < s < 1.6. This can be quan- 
tified by comparing 0.392, the fraction of spacings with 
s < 0.75, to 0.357, the corresponding fraction predicted 
using the Wigner distribution. Using the normal approx- 
imation to the binomial distribution, this discrepancy is 
nearly 7a and is thus statistically very significant (simi- 
lar conclusions are reached by the standard Kolmogorov- 
Smirnov test for comparing distributions). We conjecture 
that, as with mode intensities discussed above, the dis- 
crepancy is another manifestation of slow convergence to 
the semiclassical limit. 

One difference between our mushroom and that of Di- 
etz et. al. is that our foot supports BB orbits and theirs 
does not. Therefore to eliminate this as a cause of dif- 
ference, in Fig. [T2T i we show the ergodic NNDS with BB 
modes removed. Here BB modes were identified as those 
with fj > 0.7 but small integral on the base of the foot, 

1 /2 

namely J ' (n • r)\d n (j)j\ 2 dq < 0.1; the BB subset com- 
prises only 0.8% of the total. The difference between 
panels c) and d) is barely perceptible, indicating that BB 
modes are not a significant contribution in our setting. 

Finally, in Fig. [T3] we show the amplitude spectrum 
PQ) : — 2j=i elkjl °f the density of states, which high- 
lights contributions from periodic orbits of length 
Panel a) shows all levels, while b) and c) shows the con- 
tribution only of levels categorized as either regular or 
ergodic, according to the above method. The periodic 
peaks at the integers in panel a) (and absent in b) are 
due to the BB mode in the foot. As expected, b) con- 
tains only the regular clusters of peaks associated with 
hat orbits which unfold to polygons in the disc. Each 
cluster of peaks has an upper limit point at multiples of 
ttR = 3tt/2 corresponding to whispering-gallery rays. It 
is interesting that c) contains contributions not only from 
UPOs but from all the peaks of b) too. 



VII. CONCLUSION 

We have presented the first known high-lying eigen- 
mode calculations of Bunimovich's mushroom, which has 
unusually simple divided phase space without KAM hier- 
archy. Using a basis set adapted to the re-entrant corner, 
the Method of Particular Solutions achieves very high ac- 
curacy for low modes, and the scaling method enables us 
to find high modes orders of magnitude more efficiently 
than any other known numerical approach, allowing the 
lowest n = 16061 odd modes to be computed in rea- 
sonable time. Since statistical estimation errors scale as 
1/y/n, we are therefore able to reach the 1% level for 
many quantities. 

Chaotic modes and Husimi functions have been shown 
to be nonuniform and scarred even at mode number 
~ 45000, evidence that the semiclassical limit is reached 
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very slowly. Using a separation into regular vs chaotic 
modes, Percival's conjecture has been verified to within 
2%. A new model for dynamical tunneling (similar to 
that of Backer et. al. [Hj]) has been described, and shown 
to predict the chaotic component of predominantly- 
regular modes to within a factor of roughly an or- 
der of magnitude (over a range of 10 s ). Its predic- 
tion (via Bessel asymptotics) that the density of occur- 
rence of modes which are regular-chaotic superpositions 
dies asymptotically like E -1 / 3, agrees well with the first 
known measurement of this density. 

Our study of nearest-neighbor eigenvalue spacing finds 
good agreement with the Berry-Robnik distribution, and 
for the regular subset, good agreement with the Pois- 
son distribution. The ergodic subset shows statistically- 
significant deviations from Wigner's GOE approxima- 
tion, favoring small spacings. However we find no evi- 
dence for the dip reported at s = 0.7 by Dietz et. al. [|[; 
recall we study over 16 times their number of modes. 

This study is preliminary, and raises many interesting 
questions: Can our model for dynamical tunneling be re- 
fined to give agreement at the impressive level found in 
quantum maps 34]? Does the ergodic level-spacing dis- 
tribution eventually tend to the GOE expectation? Fi- 
nally, can spectral manifestations of stickiness [H, EH be 
detected? 
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tion || («^)""0o|l2 = n\,\/n £ N, which can be proved by in- 
duction. The Bargmann representation [52I.I53I] of a func- 
tion v : R — > C is then (ip z , v) ; the Husimi representation 
is its squared magnitude H v ^(z) := \ (tfj z , v) | 2 . We need 

a more explicit form than (|A2|) . ip z — e za ~ z a ipo fol- 
lows by the Baker-Campbell-Hausdorff formula e A+B = 
e -[A,B]/2 e A e B for [[A,B],A] = [[A,B],B] = 0. Applying 
this formula again and writing z := (qo/a + iako)/V^ 
where qo,kg el gives 



^( g )= e iWV fc °Vo(9-<Zo). 



(A3) 



This shows that the coherent state is localized in posi- 
tion (around qo) and wavenumber (around ko), thus the 
Husimi is a microlocal (phase space) represention, 



H v ,cr(qo,ko) 



v{q)e lk ^Mq-q )dq 



(A4) 



This also known as the Gabor transform or spectrogram 
(windowed Fourier transform), and it can be proven equal 
to the Wigner transform convolved by the smoothing 
function ip^. Given a normal-derivative function d n <f>j 
we periodize it in order to apply the above. We also 
scale the wavenumber by kj, thus the Birkhoff momen- 
tum coordinate is p = ko/kj. 
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APPENDIX A: HUSIMI TRANSFORM 

We define the Husimi transform [5l| of functions on K, 
for convenience reviewing the coherent state formalism 
in dimensionless (fi.-free) units. Given a width parameter 
(phase space aspect-ratio) a > 0, it is easy to show that 
the annihilation operator 

a:=-L(i + ^) (Al) 

has a kernel spanned by the L 2 -normalized Gaussian 
ip (q) := (7TCT 2 )" 1 / 4 e-« 2 / 2<T2 . We work in L 2 (R), in which 
the hermitian adjoint of a is = (q/cr — ad q )/\2. From 
the commutator [a, a'] — 1 it follows, Vz <E C, that the 
coherent state 

i, z : = e -l 2 l 2 /V°Vo (A2) 

is an eigenfunction of a with eigenvalue z. The fact that it 
is L 2 -normalized requires the Hermite-Gauss normaliza- 



APPENDIX B: INTEGRABLE PHASE-SPACE 
FRACTION 

The total phase space (restricting to the unit-speed 
momentum shell) has volume V tot = vol(0' x S 1 ) = 
27rvoliy = 27r(afe/2 + nR 2 /4). Define the function 
a(r) := 2n — 4sin~ 1 (6/2d(r)), where d(r) is the distance 
from r to the center point q c . When r is in the hat and 
d(r) £ [6/2, R], a(r) gives the measure of the set of an- 
gles in S 1 for which orbits launched from r are integrable 
(i.e. never leave the annulus rf(r) 6 [6/2, R]). The regular 
phase space volume is found by integrating a(r) over the 
quarter- annulus using polar coordinates (p,4>): 

ptt/2 nR 

U rcK = / d(p \ a(p)pdp 

JO Jb/2 

2 V 4; J b/ / 2p P 

2 / 6 6 \ 
= nR COS — Bn 

V 2R 2R J 

The same result is given without calculus using the space 
of oriented lines in a full annulus, that is, 4U rog = 2it times 
the area of the segments {(x, y) : x 2 +y 2 < R 2 , \y\ > 6/2}. 
For our parameters we get p IOS :— V IQS /V tot — 0.4549 • • • 
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We note that a related Penrose-Lifshits mushroom con- 
struction [54| continues to find use in isospectral prob- 
lems [H^] 

All calculation times are reported for one core of a 
2.4 GHz Opteron running C++ or MATLAB under 
linux/GNU 

This requirement is needed to exclude the zero-measure 
set of marginally-unstable periodic orbits (MUPOs) in 
the ergodic region which nevertheless remain in the the 
hat for all time [TJ QJ] 



FIG. 14: 20 high-eigenvalue consecutive modes, covering the 
range kj £ [499.800,499.869], with mode number j « 45000. 
Mode number increases horizontally from the top left. \4>j\ 2 
is shown with zero white and larger values darker. 




FIG. 15: Husimi distributions Hg n ^ jtCr (q,p) of the 20 high- 
eigenvalue modes shown in Fig. 1141 and in the same order. 
The q and p axes are as in Fig. |9] 



